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Abstract 

Anisotropic  surface  free  energy  is  included  in  a phase-field  model  to  determine 
the  equilibrium  shapes  of  a free  particle  in  a matrix.  Our  model  allows  for  a 
crystal-melt  interface  with  sharp  corners  due  to  a highly  anisotropic  surface 
free  energy  with  missing  orientations.  Numerical  simulations  for  various  de- 
grees of  anisotropy  were  calculated  and  they  show  excellent  agreement  with 
analytical  equilibrium  shapes. 

I.  INTRODUCTION 

Phase-field  models  provide  a numerical  technique  that  has  been  shown  to  be  useful  for  the 
study  of  solidification  processes  [1-7].  The  fundamental  component  of  a phase-field  model  is 
the  inclusion  of  an  additional  variable  <^>,  which  allows  one  set  of  thermodynamic  equations 
to  describe  a system  of  multiple  phases.  The  phase  field  is  a constant  value  in  each  bulk 
phase,  e.g.  <p  = 0 in  the  matrix  and  <p—  1 in  the  particle.  The  interface  between  phases  is 
represented  by  a smooth  transition  region  where  ef)  varies  from  zero  to  one.  Modifications 
are  made  to  the  governing  thermodynamic  functions  with  the  addition  of  gradient  energy 
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terms  that  contribute  to  the  surface  energy  of  the  system.  The  most  significant  advantage 
of  a phase-field  model  is  the  avoidance  of  explicit  tracking  of  the  interface  between  phases. 

The  incorporation  of  surface  free  energy  anisotropy  in  these  models  has  been  considered 
previously  [2,4,8-10].  The  work  of  both  Kobayashi  [2]  and  Wheeler,  Murray,  and  Schae- 
fer [9]  show  that  anisotropy  could  be  introduced  by  allowing  the  gradient  energy  coefficient  to 
depend  on  the  normal  to  the  interface.  McFadden  et  at.  [4]  followed  this  with  an  asymptotic 
analysis  for  a surface  free  energy  that  varies  smoothly  with  orientation,  i.e.  without  missing 
orientations,  and  they  recovered  the  Gibbs-Thomsonequation  in  the  sharp-interface  limit.  A 
direct  link  between  sharp  and  diffuse  surface  motion  laws  with  anisotropy  was  determined  by 
Taylor  and  Cahn  [11].  They  showed  the  dependence  of  the  energy  on  surface  normal  direction 
was  identical  with  a phase-field  or  sharp  interface  model.  Fierro  et  at.  [10]  performed  numer- 
ical simulations  of  an  anisotropic  phase-field  model  and  the  effect  which  the  anisotropy  has 
on  the  numerical  stability.  They  investigated  nonconvex  reciprocal  anisotropy  and  the  effect 
convexification  had  on  the  motion  of  a nonconserved  phase-field.  They  did  not  investigate 
a conserved  evolution  equation  or  the  resulting  equilibrium  solutions.  The  dynamics  of  the 
nonconvex  problem  should  also  be  investigated  in  correlation  with  the  numerical  technique. 

Our  aim  is  to  introduce  anisotropy  into  a simple  phase-field  model  for  a pure  crystal  in 
a liquid  melt  which  can  predict  equilibrium  shapes  with  sharp  corners.  When  the  degree 
of  anisotropy  is  high  enough,  orientations  disappear  from  the  surface,  i.e.  corners  are  de- 
veloped and  the  numerical  solution  of  the  phase-field  equations  becomes  problematic.  We 
developed  a method  that  is  numerically  stable  for  highly  anisotropic  systems  and  shows 
excellent  agreement  with  known  equilibrium  shapes. 

II.  EQUILIBRIUM  SHAPES  OF  PARTICLES  IN  A MATRIX 

The  effect  of  anisotropic  surface  free  energy  on  the  equilibrium  shape  of  two-dimensional 
domains  of  a phase  in  a matrix  can  be  described  by  the  Gibbs-Thomsonequation, 
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H = Ho  + vm( 7 + leo ) A 


where  vm  is  the  molar  volume,  A is  the  curvature,  7 is  the  surface  free  energy  and  the 
subscripts  denote  differentiation  with  respect  to  0.  the  angle  made  by  a unit  vector  normal 
to  the  interface  and  the  x axis.  Equilibrium  is  achieved  when  the  chemical  potential,  y is 
constant  throughout  the  system.  Setting  y — Hei  a constant,  equilibrium  shapes  can  be 
found  by  solving  Eq.  1 in  parametric  form  [12]; 

x = — — (7  cos  0 — 7#  sin  0) 

He  ~ Ho 

y=  — — (7  sin  0 + ye  cos  9)  (2) 

He  - Ho 

for  the  quadrant  0 < 0 < 7t/2. 

The  6 dependence  of  the  surface  free  energy  is  unspecified  to  this  point.  We  have  chosen 
the  following  7 function, 

7(0)  = 7o(l  + e4  cos  40)  (3) 

for  a crystal  with  four-fold  symmetry.  The  degree  of  anisotropy  is  set  by  the  constant 
c4  on  the  interval  [0,1).  As  the  anisotropy  is  increased,  the  crystal  shape  will  be  energy 
minimizing  when  certain  high  energy  orientations  are  missing  from  the  equilibrium  shape. 
Missing  orientations  occur  when  the  reciprocal  7 plot  first  becomes  concave  [13].  In  two 
dimensions,  concavity  of  I/7  requires  that, 

7 + lee  = 1 - 15  e4  cos 46*  = 0 (4) 

for  some  orientations.  Thus,  missing  orientations  occur  for  c4  > 1/15.  The  first,  missing 
orientation  9m , for  a given  anisotropy  can  be  determined  by  setting  y — 0 in  Eq.  2.  The 
shapes  obtained  from  Eq.  2 include  metastable  and  stable  orientations  on  the  “ears1-  that 
do  not  belong  to  the  equilibrium  crystal.  Mullins  [14]  proved  that  in  2-D  the  equilibrium 
shape  is  given  by  the  convex  shape  remaining  after  removal  of  these  “ears".  The  resulting 
equilibrium  shapes  for  various  values  of  are  shown  in  Fig.  1. 
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FIGURES 


FIG. 


(e)  (4  = 0.20 


(f)  64  = 0.50 


Equilibrium  shapes  for  anisotropic  crystals  with  various  values  of  e4. 


III.  DEVELOPMENT  OF  THE  MODEL 

Phase-field  models  are  natural  extensions  of  the  diffuse-interface  models  of  Cahn  and 
Allen  [15, 16],  Ginzburg  and  Landau  [?],  and  of  Cahn  and  Hilliard  [17, 18].  For  a single  com- 
ponent crystal,  the  phase-field  equations  are  developed  from  the  free-energy  functional  [1], 

T = fa  (fW  + j |V^|2  dCl  (5) 

on  the  region  fl.  The  energy  density  function, 

H4>)  = ~<pf  (6) 

is  a double-well  which  has  minima  at  = 0 and  4>=\.  With  this  choice  of  the  free-energy 
density,  the  barrier  height  of  the  double  well  potential  is  IT/64.  Requiring  the  free-energy 
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functional  to  decreases  monotonically  in  time  for  a conserved  variable  results  in  the  evolution 
equation, 

Dcj)  ( 8F\ 

^ = v.^/wv_j  (?) 

where  M(cf>)  is  the  mobility  and  the  ^ dependence  allows  for  a mobility  that  is  different  in 
each  phase  or  at  the  interface. 

There  are  threee  steady  state  one-dimensional  solutions  of  Eq.  7.  The  first  two  are  of 
a single-phase  a (</>=0  matrix)  or  (3  particle)  everywhere  in  the  domain.  The  third 

solution  is  for  a planar  interface, 


cp(x)  = 


1 
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1 — tanh 


(8) 


where  the  interfacial  thickness, 


(9) 


is  a balance  between  two  opposing  effects.  The  interface  tends  to  be  sharp  in  order  to 
minimize  the  regions  where  the  free-energy  density  is  positive  which  occurs  when  cp  is  between 
0 and  1.  Conversely,  the  interface  tends  to  be  diffuse  to  reduce  the  energy  associated  with 
the  gradient  of  cf).  The  phase  field  varies  from  0.1  to  0.9  over  a distance  of  AS.  The  excess 
free  energy  of  the  interface  is  related  to  the  gradient  energy  coefficient, 
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(10) 


for  the  free-energy  functional  in  Eq.  5 with  Eq.  6 [16]. 

To  describe  anisotropic  surface  energies,  we  allow  the  gradient  energy  coefficient,  e to 
depend  on  the  angle  of  the  normal  to  the  contours  of  constant  (f>  [4,5].  This  angle, 

<»> 

is  the  same  as  the  orientation  angle  described  in  Section  II.  Any  anisotropy  present  in  t will 
appear  in  the  surface  energy,  7 in  Eq.  10,  as  well  as  the  interface  thickness,  S in  Eq.  9.  We 
choose  the  anisotropy  in  e(9)  to  coincide  with  the  surface  free  energy  in  Eq.  3.  Thus, 
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e{9)  = e0(l  + e4  cos  40) 


(12) 


For  values  of  e4  > 1/15,  the  polar  plot  of  1/e  is  non-convex  when  e + egg  < 0.  This  leads 
to  an  ill-posed  evolution  equation  that  is  backward  parabolic  for  0 in  the  range  of  missing 
orientations.  The  range  of  missing  orientations  is  given  by  a common  tangent  construction 
to  the  1/e  plot,  see  Fig.  2.  Considering  the  vertical  tangent  on  the  right  side  of  the  figure, 
the  angles  at  the  tangent  points  are  extrema  in  the  abscissa  coordinate, 

d ( cos  0 


= 0 


(13) 


dO  \ e(0) 

which  reproduces  the  torque-balance  condition.  The  corner  angle  of  the  equilibrium  shape 
or  first  missing  orientation,  0m  follows  from  Eq.  13  and  therefore  stisfies, 


c( 9m ) sin  9m  "F  e#( 0m ) cos  9m  — 0 


(14) 


Performing  this  common  tangent  construction  dictates  equilibrium  or  the  balance  of  chemical 
forces  at  the  corner.  Our  method  also  removes  metastable  orientations  which  have  t + too  > 0 
but  are  not  on  the  equilibrium  shape.  This  is  a requirement  for  local  equilibrium. 


FIG.  2.  Convexifying  the  polar  plot  of  1/e  (e4  = 0.20). 


We  propose  a regularized  gradient  energy  coefficient,  e using  the  common  tangent  to 
convexify  the  1/e  plot.  The  inversion  of  the  tangent  in  Fig.  2 with  abscissa  cos(0m)  / e(6m) 
goes  into  a circle  through  the  origin  in  the  e-plane,  see  Fig.  3.  The  convexified  portion  of 
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e is  thus  a portion  of  the  circle  with  diameter  e(0m)/  cos(0m).  Restricting  attention  to  the 
quadrant  \0\  < 7t/4,  the  modification  to  e occurs  over  the  range  |#|  < 0m . Our  regularized 
gradient  energy  coefficient, 


e = 


e(0)  for  6m  < \0\  < 7r /4 

t(0m)  cos  0 


15; 


cos  (07 


for  |0|  < 6r 


will  provide  a surface  free  energy  anisotropy  that  has  the  same  equilibrium  shapes  described 
in  Section  II.  and  in  which  e + tgg  > 0.  Similarly,  the  expressions  for  the  remaining  three 
quadrants  can  be  determined,  reflecting  the  four-fold  symmetry  and  the  piece-wise  definition 
of  l. 


Using  this  convexified  gradient  energy,  Eq.  15  in  the  free-energy  functional,  Eq.  5 and 
using  Eq.  7 yields, 


ST  8f  / \ d ( _ dt  d(p\  3 (~dld<T 


(16) 


Scf)  0(p  ' ' dx  \ d9  dy  ) dy  \ dd  dx ) 

this  is  a generalized  chemical  potential.  For  non-missing  orientations  ( 9m  < |0|  < tt/4) 
Eq.  16  can  be  reformulated  into  a more  suitable  form  for  computations  [5], 


c st~  r 

j-  = ^7  ~ e2V2</i >-  et  [sin  (29)  (4>yy  - (j)xx ) + 2 cos  (20)  (j>xy] 

+ - (e/2  4-  ee")  [2  sin  (20)  4>xy  - V2</>  - cos  (26)  (<f)yy  — <pxx) 


(17) 


where  e'  = de/dO  and  e"  = d2c/d02.  For  missing  orientations  (|0|  < 6m)  Eq.  16  can  be 
simplified  greatly  due  to  the  specific  choice  of  the  convexified  gradient  energy, 


ST  df 


e(6r 


(18) 


S(f>  d<f>  \COS  {9  m)) 

for  the  quadrant  |0|  < 7r/4.  Similar  equations  can  be  derived  for  the  remaining  three  quad- 
rants. 


IV.  NUMERICAL  SIMULATIONS 

Numerical  simulations  of  the  anisotropic  particle-matrix  surface  free  energy  were  then 
performed  using  an  explicit  finite  difference  technique.  We  solve  the  evolution  equation, 
Eq.  7 on  a two-dimensional  uniform  grid  using  a finite  difference  approximation  that  is 
second  order  accurate  in  space  and  first  order  accurate  in  time.  The  mesh  spacing  is  h,  and 
the  time  step  At. 


A.  Equilibrium  Shapes 

When  the  degree  of  anisotropy  is  low  (e4  < 1/15),  the  equilibrium  shape  is  smooth  without 
corners  and  the  numerical  simulation  can  be  performed  without  difficulty.  For  the  spatial 
derivatives,  simple  centered  finite  differencing  formula  can  be  used  throughout  the  domain. 
The  equilibrium  shapes  that  were  calculated  numerically  for  t4  = 0.06  match  well  with  the 
analytical  shapes  calculated  by  Eq.  2 in  Section  II.  Fig.  4 shows  this  agreement  in  both 
the  numerical  and  the  analytical  equilibrium  shapes  (shown  by  the  dashed  and  solid  lines 
respectively).  In  Fig.  4(a)  the  dashed  numerical  equilibrium  shape  is  constructed  from  a 
contour  of  constant  phase-field  parameter,  (f)  = 0.5.  In  Fig.  4(b)  is  a polar  plot  of  the  surface 
normal  angle  as  a function  of  the  polar  angle  of  the  point  on  the  interface  at  which  the 
normal  angle  is  measured. 
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(a)  equilibrium  shapes  (b)  surface  normal  orientations 

FIG.  4.  Comparison  of  numerical  dashed  line  and  analytical  solid  line  (a)  equilibrium  shapes 
and  (b)  surface  normal  orientations,  6 for  particles  with  non-missing  surface  energy  anisotropy 
(€4  = 0.06). 


For  anisotropies  with  e4  > 1/15,  orientations  disappear  from  the  equilibrium  shape  and 
the  numerical  calculation  becomes  more  involved.  First,  special  care  is  taken  to  accurately 
evaluate  0,  the  angle  of  the  normal  to  the  contours  of  constant  <f>  given  by  Eq.  11.  The 
first  derivatives  of  <p  across  corners  are  erroneously  calculated  by  centered  finite  differencing 
formula.  For  instance,  if  a mesh  point  is  located  exactly  on  a corner,  then  0 = 0 using  a 
centered  formula  (he.  Eq.  19  and  Fig.  5(a)).  However,  a one-sided  formula  (he.  Eq.  20  and 
Fig.  5(b))  would  yield  the  correct  angle  on  either  side  of  the  corner.  Although  this  difficulty 
can  be  avoided  by  rotating  the  mesh  grid  relative  to  the  surface  anisotropy  by  45°  which 
would  reduce  the  inaccuracy  of  the  centered  differencing,  this  does  not  completely  solve  the 
problem  especially  for  very  high  anisotropies. 


j.  + i)  - 4>{iJ  - 1) 



-3 <j){ij)  + 4 <f)(i,j  + 1)  - 0(i,  j + 2) 

Wffi)  - 777 


(19) 

(20) 
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(a)  centered  (b)  one-sided 

FIG.  5.  Diagram  of  6 calculation  by  the  related  first  derivatives  using  (a)  centered  and  (b) 
one-sided  differencing  formula  across  the  corner,  (where  0C  is  the  corner  angle)  IA  the  x direction, 
centered  differencing  is  used. 

The  second  derivatives  (^>xx,  (j)yy,  (foxy)  are  also  inaccurate  using  centered  differencing  for- 
mulae near  a corner.  Unfortunately,  using  a similar  approach  as  that  used  in  evaluating 
the  first  derivatives  yields  a method  that  is  unstable  in  time.  Ignoring  any  error  caused  by 
using  centered  second  derivatives,  we  are  still  able  to  calculate  very  good  equilibrium  shapes. 
Shown  in  Fig.  6 is  a comparison  of  the  analytical  and  numerical  crystal  shapes  for  e4=0.08 
(which  is  greater  than  1/15  and  thus  has  corners). 
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(a)  equilibrium  shapes  (b)  surface  normal  orientations 

FIG.  6.  Comparison  of  numerical  dashed  line  and  analytical  solid  line  (a)  equilibrium  shapes  and 
(b)  surface  normal  orientations  (6)  for  particles  with  missing  surface  energy  anisotropy  (c4  = 0.08). 

Remarkably,  these  numerical  shapes  have  discontinuous  chemical  potentials  at  the  cor- 
ners, but  are  still  numerically  stable.  Using  centered  second  derivatives  results  in  a corner 
chemical  potential  that  changes  from  small  to  large  values  nearly  every  iteration.  Completely 
one-sided  second  derivatives  are  accurate,  but  unstable.  Accurately  resolving  the  chemical 
potential  while  also  maintaining  the  stability  of  the  numerical  method  is  the  trade-off  that 
must  be  balanced.  Thus,  a mixed  differencing  scheme  was  developed  for  more  accurate  cal- 
culations near  these  sharp  corners.  A comparison  between  the  chemical  potentials  calculated 
with  this  mixed  differencing  scheme  verses  that  of  a simple  centered  one  is  shown  in  Fig.  7. 
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FIG.  7.  Polar  polat  of  the  chemical  potential  along 
and  (b)  mixed  second  derivative  differencing  formula. 


(b)  mixed 

the  interface  for  €4  = 0.08  using  (a)  centered 


In  our  mixed  differencing  scheme,  the  first  issue  is  to  locate  the  corners  in  our  phase-field 
based  on  the  value  of  0 and  then  choose  the  differencing  scheme  for  the  second  derivatives 
that  go  across  each  corner.  If  \0\  <0m , then  the  convexified  gradient  energy  is  used  (Eq.  15) 
and  the  chemical  potential  is  determined  by  Eq.  18.  The  chemical  potential  only  depends 
on  the  second  x derivative  at  the  corner,  but  not  across  it;  hence,  centered  differencing  is 
appropriate.  For  |0|£0m,  the  full  equation  for  the  chemical  potential  (Eq.  17)  must  be 
used.  In  this  region,  a weighted  average  of  one-sided  and  slanted  differencing  (away  from  the 
corner)  is  used  for  those  derivatives  across  the  corner.  The  weighted  average  is  95  to  99.9% 
one-sided  depending  on  the  anisotropy.  If  0m  < \6\  < 7r/8  -f-  9ml 2,  then  a slanted  formula  is 
used.  Finally  for  7t/8  + 0m/2  < |#|  < 7t/4,  centered  differencing  is  used  for  all  derivatives 
because  it  is  significantly  far  away  from  the  sharp  corner.  Examples  of  differencing  formula 
for  centered,  slanted,  and  one-sided  are  given  in  Eqs.  21-23. 


^yy  ( b J ) 


(t)yy(li  J ) 


~ 1)  - + 1) 
h2 

11  - 1)  - 20 <f>(i,j)  + + 1)  + 4 + 2)  - + 3) 

12  h2 


(21) 


221 
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(23) 


. 2 <f>(i,j)  - 5 <f>(ij  + 1)  + 4 <f>(i,j  + 2)  - + 3) 

(Pyy\‘l’>  J ) — 1^2 

A similar  algorithm,  based  on  the  value  of  0,  for  the  remaining  three  quadrants  is  straight- 
forward. 

The  chemical  potential  computed  with  the  centered  differencing  algorithm  is  discontin- 
uous but  stable,  due  to  the  regularization  of  the  gradient  energy  coefficient.  Switching  to  a 
mixed  differencing  algorithm  results  in  a smooth  constant  chemical  potential  (Fig.  7).  How- 
ever, the  equilibrium  shape  is  unchanged.  Computing  good  equilibrium  shapes  for  values 
of  e4  < 0.5  is  possible  with  either  differencing  technique.  This  limit  in  surface  free  energy 
anisotropy  is  not  due  to  the  finite  difference  approximation,  but  instead  a result  of  the  mesh 
resolution.  Shown  in  Figs.  8 and  9 compare  the  analytical  and  numerical  crystal  shapes  for 
64  = 0.15  and  £4  = 0.30,  respectively. 


(a)  equilibrium  shapes 


(b)  surface  normal  orientations 


FIG.  8.  Comparison  of  numerical  dashed  line  and  analytical  solid  line  (a)  equilibrium  shapes  and 
(b)  surface  normal  orientations  ( 6 ) for  particles  with  missing  surface  energy  anisotropy  (£4  = 0.15). 
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(a)  equilibrium  shapes  (b)  surface  normal  orientations 

FIG.  9.  Comparison  of  numerical  dashed  line  and  analytical  solid  line  (a)  equilibrium  shapes  and 
(b)  surface  normal  orientations  (0)  for  particles  with  missing  surface  energy  anisotropy  — 0.30). 

B.  Dynamics 

The  dynamics  of  forming  crystal  shapes  with  corners  from  an  initially  smooth  surface 
was  investigated  with  a centered  differencing  scheme.  The  mixed  scheme  presented  above 
is  optimized  for  shapes  that  have  corners  and  is  inadequate  for  a circular  initial  condition. 
The  regions  on  the  initial  shape  in  the  missing  range  of  orientations  evolve  slower  than  those 
just  outside  of  this  area.  The  evolution  is  shown  in  Fig.  10  by  surface  contours  (<f>  = 0.5)  at 
various  times  shown  atop  each  other.  This  phenomenon  is  non-physical  and  is  most  likely 
a result  of  the  regular  zed  gradient  energy.  However,  Fig.  10  illustrates  that  the  method  is 
stable  even  though  the  initial  shape  had  orientations  for  which  e + tee  < 0. 
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FIG.  10.  Surface  contours  during  computation  of  equilibrium  shape  (e  = 0.15). 

Once  the  corner  has  formed,  its  motion  is  not  affected  by  the  differencing  scheme.  This 
was  tested  by  tracking  the  corner  positions  during  the  evolution  from  a rectangular  shape  to 
a shape  with  equal  length  sides  (Fig.  11).  No  significant  difference  was  found  between  the 
computations  using  centered  and  mixed  differencing  formula.  Corner  motion  is  unaltered  by 
the  discontinuity  of  the  chemical  potential  present  with  the  centered  scheme. 


FIG.  11.  Diagram  of  corner  motion  test,  showing  the  initial  and  final  surface  shapes  (c  = 0.15). 

Any  difference  between  the  shapes  given  by  a centered  or  mixed  differencing  scheme  is 
restricted  to  one  or  two  mesh  points  around  each  corner.  These  slight  differences  are  not 
apparent  in  the  equilibrium  shapes  shown,  due  to  the  number  of  mesh  points  used  for  each 
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computation  (mesh  261  x 261).  With  a centered  calculation,  the  corners  are  rounded  by 
orientations  that  should  not  appear  on  the  equilibrium  shape.  However,  there  appears  to 
be  no  other  significant  disadvantage  with  using  simple  centered  differencing  for  the  second 
derivatives;  the  time  step  stability  is  the  same  for  either  method. 

A computation  was  performed  on  a sinusoidal  interface  to  demonstrate  the  flexibility  of 
this  phase-field  method,  see  Fig  12.  The  initially  smooth  surface  was  imposed  with  a surface 
energy  anisotropy  resulting  in  corners.  Initially  the  jumps  in  normal  angle  or  corners  develop 
very  quikly.  Along  the  steep  slides  new  corners  appear  (t  — b)  due  to  the  enforcement  of 
local  equilibrium  required  by  our  regularization  of  the  gradient  energy.  These  corners  then 
coarsen  (t  = c)  and  disapear  (/  = d)  as  the  phase-field  evolves.  The  final  interface  only 
contains  stable  orientations.  However,  it  is  metastable  because  the  curved  surface  area  can 
be  reduced  by  coarsening.  This  computation  shows  that  the  method  can  both  introduce  new 
and  remove  old  corners  without  any  implicit  tracking  of  the  interface. 


FIG.  12.  Surface  contours  {4>  — 0.5)  at  various  times  for  a periodic  sinusoidal  computation  for 
c4  = 0.15. 
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V.  CONCLUSION 


We  have  introduced  anisotropic  surface  free  energy  into  a phase-field  model  and  com- 
puted the  equilibrium  shapes  for  crystals  both  with  and  without  sharp  corners.  We  have 
shown  excellent  agreement  between  our  equilibrium  shapes  and  the  analytical  ones  for  highly 
anisotropic  particles.  With  the  numerical  solution  of  the  phase-field  equations  for  this  sim- 
ple, two  dimensional,  crystal-melt  system,  we  were  able  to  determine  the  equilibrium  shapes 
for  various  degrees  of  anisotropy.  No  obvious  disadvantage  using  simple  centered  differencing 
for  the  second  derivatives  has  been  determined.  Numerical  stability,  dynamics,  and  most 
importantly  equilibrium  shape  are  unaffected  by  the  discontinuity  of  the  chemical  potential 
present  with  the  centered  algorithm.  This  validated  phase-field  method  allows  similar  for- 
mulations, for  surface  energy  anisotropy,  to  be  included  into  more  sophisticated  phase-field 
simulations. 
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